pomp (version 1.6)

Particle Markov Chain Monte Carlo: The particle Markov chain Metropolis-Hastings algorithm

Description

The Particle MCMC algorithm for estimating the parameters of a partially-observed Markov process. Running pmcmc causes a particle random-walk Metropolis-Hastings Markov chain algorithm to run for the specified number of proposals.

Usage

"pmcmc"(object, Nmcmc = 1, start, proposal, Np, tol = 1e-17, max.fail = Inf, verbose = getOption("verbose"), ...)
"pmcmc"(object, Nmcmc = 1, Np, tol, ...)
"pmcmc"(object, Nmcmc, start, proposal, Np, tol, max.fail = Inf, verbose = getOption("verbose"), ...)
"continue"(object, Nmcmc = 1, ...)

Arguments

object

An object of class pomp.

Nmcmc

The number of PMCMC iterations to perform.

start

named numeric vector; the starting guess of the parameters.

proposal

optional function that draws from the proposal distribution. Currently, the proposal distribution must be symmetric for proper inference: it is the user's responsibility to ensure that it is. Several functions that construct appropriate proposal function are provided: see MCMC proposal functions for more information.

Np

a positive integer; the number of particles to use in each filtering operation.

tol

numeric scalar; particles with log likelihood below tol are considered to be “lost”. A filtering failure occurs when, at some time point, all particles are lost.

max.fail

integer; maximum number of filtering failures permitted. If the number of failures exceeds this number, execution will terminate with an error.

verbose

logical; if TRUE, print progress reports.

...

additional arguments that override the defaults.

Value

An object of class pmcmc.

Re-running PMCMC Iterations

To re-run a sequence of PMCMC iterations, one can use the pmcmc method on a pmcmc object. By default, the same parameters used for the original PMCMC run are re-used (except for tol, max.fail, and verbose, the defaults of which are shown above). If one does specify additional arguments, these will override the defaults.

Continuing PMCMC Iterations

One can continue a series of PMCMC iterations from where one left off using the continue method. A call to pmcmc to perform Nmcmc=m iterations followed by a call to continue to perform Nmcmc=n iterations will produce precisely the same effect as a single call to pmcmc to perform Nmcmc=m+n iterations. By default, all the algorithmic parameters are the same as used in the original call to pmcmc. Additional arguments will override the defaults.

Details

pmcmc implements an MCMC algorithm in which the true likelihood of the data is replaced by an unbiased estimate computed by a particle filter. This gives an asymptotically correct Bayesian procedure for parameter estimation (Andrieu and Roberts, 2009). Note that pmcmc does not make use of any parameter transformations supplied by the user.

Methods

c

Concatenates pmcmc objects into a pmcmcList.

conv.rec(object, pars)

returns the columns of the convergence-record matrix corresponding to the names in pars as an object of class mcmc or mcmc.list.

filter.traj(object, vars)

returns filter trajectories from a pmcmc or pmcmcList object.

plot

Diagnostic plots.

logLik

Returns the value in the loglik slot.

coef

Returns the last state of the MCMC chain. As such, it's not very useful for inference.

covmat(object, start, thin, expand)

computes the empirical covariance matrix of the MCMC samples beginning with iteration start and thinning by factor thin. It expands this by a factor expand^2/n, where n is the number of parameters estimated. By default, expand=2.38. The intention is that the resulting matrix is a suitable input to the proposal function mvn.rw.

References

C. Andrieu, A. Doucet and R. Holenstein, Particle Markov chain Monte Carlo methods, J. R. Stat. Soc. B, to appear, 2010.

C. Andrieu and G.O. Roberts, The pseudo-marginal approach for efficient computation, Ann. Stat. 37:697-725, 2009.

See Also

pomp, pfilter, MCMC proposal distributions, and the tutorials on the package website.

Examples

## Not run: 
# library(pomp)
# 
# pompExample(ou2)
# 
# pmcmc(
#       pomp(ou2,dprior=Csnippet("
#    lik = dnorm(alpha_2,-0.5,1,1) + dnorm(alpha_3,0.3,1,1);
#    lik = (give_log) ? lik : exp(lik);"),
#            paramnames=c("alpha.2","alpha.3")),
#       Nmcmc=2000,Np=500,verbose=TRUE,
#       proposal=mvn.rw.adaptive(rw.sd=c(alpha.2=0.01,alpha.3=0.01),
#         scale.start=200,shape.start=100)) -> chain
# continue(chain,Nmcmc=2000,proposal=mvn.rw(covmat(chain))) -> chain
# plot(chain)
# chain <- pmcmc(chain)
# plot(chain)
# 
# library(coda)
# trace <- window(conv.rec(chain,c("alpha.2","alpha.3")),start=2000)
# rejectionRate(trace)
# effectiveSize(trace)
# autocorr.diag(trace)
# 
# summary(trace)
# plot(trace)
# 
# heidel.diag(trace)
# geweke.diag(trace)
# ## End(Not run)